This notebook computes the naturalness areas for each point in a 'sample_points' vector map based on a 'naturalness' vector map. The computed areas are exported to a csv file.
This notebook uses GRASS GIS (7.0.4), and must be run inside of a GRASS environment (start the jupyter notebook server from the GRASS command line).
The code in this notebook was modified from examples in "How to write a Python GRASS GIS 7 addon" developed for a FOSS4G Europe 2015 workshop. The code is contained in the python-grass-addon repository on GitHub.com.
naturalness – vector map representing naturalness
naturalness_value_field – name of the field containing the naturalness values
sample_points – vector map containing observer points
radius – radius in map units around each sample point to calculate land cover areas within
use_viewshed – consider only 'visible' land cover, based on a viewshed from each sample point
In [ ]:
naturalness = 'naturalness'
In [ ]:
naturalness_value_field = 'value'
In [ ]:
n_types = 7
In [ ]:
sample_points = 'sample_points_field'
In [ ]:
radius = 500
In [ ]:
use_viewshed = False
get viewshed suffix (for filenames)
In [ ]:
def getViewshedSuffix():
if use_viewshed:
viewshed_suffix = '_viewshed'
else:
viewshed_suffix = ''
return viewshed_suffix
In [ ]:
naturalness_area_table_filename = "/home/ubuntu/naturalness_areas_{0}m{1}.csv".format(radius, getViewshedSuffix())
In [ ]:
import numpy as np
import pandas
import pyprind
GRASS import statements
In [ ]:
import grass.script as gscript
In [ ]:
from grass.pygrass.vector.geometry import Point
from grass.pygrass.vector import Vector
from grass.pygrass.vector import VectorTopo
from grass.pygrass.vector.table import DBlinks
connect to attribute table
In [ ]:
def connectToAttributeTable(map):
vector = VectorTopo(map)
vector.open(mode='r')
dblinks = DBlinks(vector.c_mapinfo)
link = dblinks[0]
return link.table()
extract point from vector map
In [ ]:
def extractPoint(input, ID, output):
where = 'ID = {0}'.format(ID)
type = 'point'
gscript.read_command('v.extract',
input=input,
where=where,
output=output,
type=type,
overwrite=True)
create buffer around point
In [ ]:
def bufferPoint(input, output, radius):
gscript.read_command('v.buffer',
input=input,
output=output,
type='point',
distance=radius,
overwrite=True)
convert viewshed from raster to vector map
In [ ]:
def vectorizeViewshed(input, ID, output):
type = 'area'
column = 'visible'
gscript.read_command('r.to.vect',
input=input,
output=output,
type=type,
column=column,
overwrite=True)
overlay a vector map on an underlying vector map using 'and' selection operator
In [ ]:
def overlay(overlay, underlay, output):
operator='and'
gscript.read_command('v.overlay',
ainput=overlay,
binput=underlay,
operator=operator,
output=output,
overwrite=True)
add area column to vector map
In [ ]:
def calculateAreas(map):
#add new area column
gscript.read_command('v.db.addcolumn',
map=map,
columns="area_square_meters DOUBLE PRECISION")
#compute area and insert into area column
gscript.read_command('v.to.db',
map=map,
type='centroid',
option='area',
columns='area_square_meters',
unit='meters')
create table showing total area by landcover type
In [ ]:
def getNaturalnessAreaByValue(map):
#get area data
table = connectToAttributeTable(map=map)
table.filters.select()
columns = table.columns.names()
cursor = table.execute()
result = np.array(cursor.fetchall())
cursor.close()
data = pandas.DataFrame(result, columns=columns).set_index('cat')
#make sure naturalness_value_field is a numeric data type
data['b_' + naturalness_value_field] = pandas.to_numeric(data['b_' + naturalness_value_field])
#create naturalness column
data['naturalness'] = np.nan
#define naturalness value categories (round values)
for index, row in data.iterrows():
naturalness = np.round(row['b_' + naturalness_value_field])
data.set_value(index, 'naturalness', int(naturalness))
#sum areas by naturalness value
data['area_square_meters'] = pandas.to_numeric(data['area_square_meters'])
areas = data[['naturalness', 'area_square_meters']].groupby(by='naturalness').sum()
#calculate mean
total_area = data['area_square_meters'].sum()
percent_area = data['area_square_meters'] / total_area
weighted_area = percent_area * data['b_' + naturalness_value_field]
mean = weighted_area.sum()
#add to areas dataframe
areas = areas.set_value('mean', 'area_square_meters', mean)
return areas
export vector map to a shapefile
In [ ]:
def exportVectorToShapefile(map, output):
gscript.read_command('v.out.ogr',
input=map,
format='ESRI_Shapefile',
output=output,
flags='e',
overwrite=True)
get info about a vector map
In [ ]:
def getVectorMapInfo(map):
return gscript.read_command('v.info', map=map)
connect to 'sample_points' attribute table
In [ ]:
point_table = connectToAttributeTable(map=sample_points)
point_table.filters.select()
columns = point_table.columns.names()
cursor = point_table.execute()
result = np.array(cursor.fetchall())
cursor.close()
point_data = pandas.DataFrame(result, columns=columns).set_index('cat')
loop through sample points
In [ ]:
with Vector(sample_points, mode='r') as points:
#setup progress bar
progress_bar = pyprind.ProgBar(points.n_lines, bar_char='█', title='Naturalness analysis progress', monitor=True, stream=1, width=50)
#iterate through points
for point in points:
#get point ID (SiteID)
ID = point_data['ID'][point.cat-1]
#update progress bar
progress_bar.update(item_id=ID)
#buffer current point
extractPoint(input='sample_points_field', ID=ID, output='tmp_buffer_point')
bufferPoint(input='tmp_buffer_point', output='tmp_point_buffer', radius=radius)
#set buffer as overlay
overlay_input = 'tmp_point_buffer'
#consider only visible naturalness if 'use_viewshed' = True
if use_viewshed:
viewshed = 'vect_{0}_viewshed'.format(ID)
visible_viewshed = 'vect_{0}_viewshed_{1}m'.format(ID, radius)
#vectorize viewshed
vectorizeViewshed(input='{0}_viewshed'.format(ID), ID=ID, output=viewshed)
#overlay buffer on viewshed
overlay(overlay=overlay_input,
underlay=viewshed,
output=visible_viewshed)
#set overlay to the visible viewshed
overlay_input = visible_viewshed
overlay_output = 'vect_{0}_naturalness_{1}m{2}'.format(ID, radius, getViewshedSuffix())
#overlay naturalness
overlay(overlay=overlay_input,
underlay=naturalness,
output=overlay_output)
#calculate naturalness area
calculateAreas(map=overlay_output)
create table with naturalness areas by type for each point
In [ ]:
#create table
index_start = 0
''' set the first index number,
allowing easier insertion into a database table that already contains
area calculations with other parameters
(i.e radius and use_viewshed)'''
columns = ['ID', 'SiteID', 'IncludedArea']
columns = columns + [ str(n) for n in range(1,n_types+1) ] + ['mean']
area_table = pandas.DataFrame(columns=columns)
#set naming variables
included_area = '{0}m{1}'.format(radius, getViewshedSuffix()) #0=radius, 1=viewshed_suffix
map_pattern = 'vect_{0}_naturalness_{1}m{2}' #0=ID, 1=radius, 2=viewshed_suffix
#iterate through points
for index, point in point_data.iterrows():
ID = point['ID']
map = map_pattern.format(ID, radius, getViewshedSuffix())
#initiate row
row = {'ID':"{0:.3g}".format(int(index) + index_start),
'SiteID': str(ID),
'IncludedArea': included_area}
#get naturalness areas
areas = getNaturalnessAreaByValue(map)
#iterate through area types
for index, area in areas.iterrows():
#add area to row
try:
row["{0:.3g}".format(int(index))] = area['area_square_meters']
except ValueError:
row[index] = area['area_square_meters']
#append row to table
area_table = area_table.append(row, ignore_index=True)
area_table.set_index('ID', inplace=True)
#export table to file
area_table.to_csv(naturalness_area_table_filename, header=False)